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Abstract 

Drastic change in dynamics and statistics in a chemical reaction system, induced by smallness 
in the molecule number, is reported. Through stochastic simulations for random catalytic reaction 
networks, transition to a novel state is observed with the decrease in the total molecule number 
N, characterized by: i) large fluctuations in chemical concentrations as a result of intermittent 
switching over several states with extinction of some molecule species and ii) strong deviation 
of time averaged distribution of chemical concentrations from that expected in the continuum 
limit, i.e., N — > oo. The origin of transition is explained by the deficiency of molecule leading to 
termination of some reactions. The critical number of molecules for the transition is obtained as a 
function of the number of molecules species M and that of reaction paths K, while total reaction 
rates, scaled properly, are shown to follow a universal form as a function of NK/M. 

PACS numbers: 87.16.Yc, 82.39.Rt, 05.40.-a 
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1. INTRODUCTION 



In intra- cellular biochemical reaction processes, some chemical species often play an 
important role at extremely low concentrations, amounting to only a few molecules per 

ceiiflsy,^. In such situations, the fluctuations and discreteness in the molecule number 

are obviously important. On the other hand, in rate equations, generally adopted in chem- 
ical kinetics, the concentration of each chemical species is treated as a continuous variable, 
and the fluctuations and the discreteness of the number of molecules are neglected. However, 
if the molecule number is not very large, the number fluctuations as well as discreteness in 
the number, rather than a continuous concentration, has to be considered seriously. 

Of course, effects of fluctuations in concentrations are considered by using stochas- 
tic differential equations. Indeed, several non-trivial noise-induced phenomena have been 
reported^, [f3, Q|. However, in most of such studies, discreteness in the molecule numbers, 
i.e., the number being 0,1,2,.., has not been considered seriously. 

Recently, Togashi and Kaneko reported a drastic change in the steady distribution of 
chemical concentrations as a result of discreteness in the molecule number, by studying 
a catalytic reaction network with a few molecules species 

a a. 

Novel types of dissipa- 

tive structure formation, induced by very low concentration molecules, have also been 



investigated 
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15| in a class of reaction-diffusion systems or the models 



of biochemical reactions in cells. The observed novel states in these studies are a result 
of fluctuations and discreteness in the molecular numbers, in particular of extinction and 
re-emergence of some molecule species, which alternate in time. 

Relevance of such discreteness in molecule number to emergence of novel states should not 
be restricted to a simple reaction network with a few molecule species, but is also expected 
to exist in a wide variety of chemical reaction systems with a large number of species. 
Here, dynamics and statistics of chemical reaction systems with a finite (small) number of 
molecules and a large number of molecule species have to be investigated. Such study is 
important not only for biochemical reaction kinetics generally, but also as a problem of non- 
equilibrium statistical mechanics. Here, the criterion on "smallness" in the number itself 
should be clarified as a condition for discreteness-induced transition for a given chemical 
reaction network system. 

In this paper, we study discreteness-induced transition in a reaction network where a 
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large number of chemical species is connected by catalytic reaction pat hs chosen randomly. 
Use of random catalytic reaction network is pioneered by Kauffman [id \v\ for the problem 
of the origin of life, while studies in a growing cell model consisting of such reaction network 
have unveiled universal statistical behaviors of chemical concentrations, which are confirmed 



in the gene expression data in the present cells [lg, Il9| . Here, we are interested in how 
discreteness in the molecule numbers affect global behavior of chemical reaction dynamics 
in such network. For simplicity, we only consider the reaction network whose steady state is 
unique and stationary when the number of molecules is infinite, i.e., the corresponding rate 
equation has just a unique fixed point attractor. Even in such a simple system, we find the 
following transition to a discreteness-induced state, which appears when the total number 
of the molecules is below a critical value; 

A) Chemical concentrations exhibit intermittent switching among several states with 
distinct chemical compositions. 

B) The long time average of the chemical concentrations deviate distinctly from those 
expected in the continuum limit with a large number of molecules where the rate equation 
description is valid. 

We will also obtain the critical molecule number for this discreteness induced transition, 
whose dependence on the number of chemical species and path ratio of the reaction network 
will be derived. 

In the next section, we introduce a specific reaction network model, while numerical 
results to show the above transition with the decrease in the molecule number are presented 
in section 3. Possible mechanism for the transition is described in section 4, where deficiency 
in some molecules is shown to introduce switching among several effective reaction networks 
that consist only of non-vanishing chemical species. In section 5, we obtain the critical 
number of molecules for the transition as a function of the number of molecule species and 
connectivity in the network. Summary and discussion will be given in section 6. 

2. MODEL 

Now, we introduce a simple model of a network of elementary reactions that consists 



of a variety of chemical species lg]. State of the system is represented by a set of num- 



bers (ni,ri2, ....,njvf), where nA= 0, 1, ...) indicates the number of molecules of the chemical 
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species i(l < i < M), with M as a total number of molecule species. Here, the total number 
of molecules is fixed at N, and accordingly < rii < N. For the chemical reaction dynam- 
ics, we choose a catalytic network among these chemical species, where each reaction from 
a chemical B to another chemical A is assumed to be catalyzed by a third chemical C; i.e., 
B + C — > A + C (see Fig. 1). The reaction coefficients are set to be identical for all reactions 
for simplicity, and chosen to be 1. Then, the growth rate in ua (or the decay rate in ub) 
through this reaction is give by riBnc/N 2 , on the average. 

The connection paths in a reaction network are chosen randomly (and then fixed), where 
the average number of the reaction paths from a chemical i to any other chemical species 
j catalyzed by a chemical I is set at a given connection number K. We do not include 
auto-catalytic reaction in the form of B + C — > 2C, because such type of reaction is not 
usually elementary but is realized as a result of a series of (non-auto-catalytic) elementary 
reactions. Also, inclusion of auto-catalytic paths sometimes leads to non-fixed-point or 
multiple attractors, which makes discussion on the discreteness effect complicated. 

We also include a flow of chemicals into and out of the system from the reservoir. With 
this process, the total number of molecules is fixed so that a molecule is regarded to be 
replaced by some other with a certain rate. Instead of considering such flow, one can 
equivalently consider a combination of decomposition and synthesis of some chemicals, or 
non-catalytic changes between different molecule species that are chosen randomly from all 
chemicals with equal probability. For simplicity, we assume that this 'non-catalytic' change 
occurs with equal probability e for all molecules, while its rate is much smaller than that of 
the catalytic reaction. 

Numerical simulations are carried out by iterating the following stochastic processes. 
First, we randomly pick up a pair of molecules and if the pair is a substrate and catalyst 
according to the reaction network, the substrate molecule is transformed to the product 
molecule according to the reaction. Second, we randomly pick up a molecule and transform 
it by a non-catalytic change, with a given, much lower, rate. Here, a unit time is given as 
the time span in which Monte-Carlo steps for catalytic reactions are repeated N/2 times and 
those for non-catalytic processes are repeated N times. In each time, each molecule collides 
with another molecule once on average, to check if the catalytic reaction occurs, while 
it is transformed transformed to some other chemical with probability e by non-catalytic 
process. Numerically, we apply this stochastic simulation, while in the cases with N » M 
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(i.e., continuous limit), the reaction dynamics are represented by the following rate equation, 

jM j',V 

where C(i,j, I) is 1 if there is a reaction i + I — > j + I, and otherwise. 

In the following sections, we present numerical results of the stochastic simulations and 
show how the steady state properties of such catalytic reaction networks depend on the 
number of chemical species M and of molecules N. In this paper, we consider the case with 
K as K c < K « M where K c indicates the critical connection number of the percolation 
transition in random networks. It is noted that the rate equation (1) in the continuum limit 
has only a unique fixed point attractor that gives the concentrations at the steady state. 

3. DISCRETENESS-INDUCED TRANSITION 

We show the results of the stochastic simulation of catalytic reaction networks with 
e = 10~ 4 . We study how dynamical aspects of the system change depending on the number 
of molecules N, by taking a given reaction network. Figures 2(a) and (b) show typical 
temporal evolutions of the concentrations of each chemical species rii for a reaction network 
with M = 100 and K = 12. Two cases with (a) large (N = 800) and (b) small (N = 12) 
iV are plotted, while Figure 2(c) shows the reaction ratio (RR) of these cases, where RR 
is defined by [Number of reacting molecules] /iV at each time. If N is much larger than M 
(e.g., N = 800), as seen in Figs. 2(a) and (c), the distribution of rij is almost stationary and 
RR takes a constant value except for small fluctuations. 

On the other hand, if N is much smaller than some value (that is of the order of M) 
(e.g., N = 12), the distribution of rii exhibits remarkable changes in time, and RR fluctuates 
intermittently between large values and 0, as seen in Figs. 2(b) and (c). Such non- stationary 
behavior is caused by the discreteness in the molecular number. When N is much smaller 
than some threshold of the order of M, the concentrations of some chemical species n» go 
to at some time instance. Then the number of reactive chemicals or of the catalysts often 
goes to 0, for all reaction paths. Then, the catalytic reactions of the system freeze, while 
the system can escape from such freezing state as a result of non-catalytic changes (flow of 
molecules). Thus, RR changes intermittently with time. 

Next, we focus on typical statistical aspects of the system by changing N. Long time 
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average of the distributions of the chemical concentrations < rii/N > for a typical reaction 
network with M = 100 and K = 12 is plotted in Fig. 3(a) (for large N, N = 800) and 
Fig. 3(b) (for a small N, N = 12). As seen in these figures, the profiles of the average 
distribution are quite different between the two. In Figs. 3(c) and (d), we show < rii/N > 
and the rank of abundances of each chemical concentration, which is labeled in the order of 
magnitude of < rii/N >. These are plotted as a function of N for a given reaction network, 
where each successive curve indicates the change in < rii/N > and the rank of each chemical 
concentration. 

The results in Figs. 3(c) and (d) suggest the existence of a critical value of N, denoted by 
N c . For N smaller than iV c , the chemical abundances of each species or their rank changes 
sensitively with N. The variance of < rii/N > over time takes a maximum at some N slightly 
smaller than this critical value N c . On the other hand, for N larger than N c , < rii/N > 
are almost constant except for the small fluctuations, and indeed the profile of < rii/N > 
for such larger N agrees with that obtained by the rate equation (1), i.e., the system is well 
described by the continuum limit N — > oo. 

Now, we focus on the transition of reaction dynamics at N c . For N larger than N c , the 
temporal evolution of the distribution of is almost stationary except for small fluctuations 
as in Fig. 2(a). On the other hand, for iV < N c , there appears intermittent switching between 
different distributions of as shown in Fig. 2(b). In this case, the total reaction ratio is 
much smaller and the temporal fluctuation of each chemical concentration is much larger 
than the case with iV > N c . We have computed N dependence of the average reaction ratio 

< RR > and the average temporal fluctuation (ATF) of all chemical concentrations. Here, 

< RR > is defined as the long time average of RR, and ATF is defined as the long time 
average of jjJ2i( n i— < n i >) 2 - 

Figure 3(e) shows < RR > and ATF/(JV/M), plotted as a function of N for the same 
reaction network. In this figure, < RR > starts to decrease drastically with the decrease in 
N, at the above mentioned critical value. This figure also shows that ATF is proportional to 
N for N > N c , while ATF /(JV/M) increases sharply with the decrease in N for N < N c . This 
increase in the fluctuation is consistent with the fact that the switching behavior becomes 
dominant for N < N c . 

As shown in Fig. 1, switching over several states with different effective network occurs 
for small N, as will be also discussed in the next section. As N is decreased further below 
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N c , the reaction occurs rarely, and thus the frequency of such switching decreases, while 
for large N, the dynamics exhibit only small fluctuations around a stationary state. Thus, 
the fluctuation is expected to have a peak around N = N c where the switching occurs most 
frequently. 

4. EFFECTS OF MOLECULAR DEFICIENCY IN A SMALL REACTION NET- 
WORK 

The behavior observed in the previous section is rather common over a variety of catalytic 
reaction networks. The dynamics and the time averaged distributions of rii/N for a case 
with iV << N c differ distinctly from those with iV — > oo, because of the discreteness in the 
number of molecules. Here we discuss the mechanism of this transition, by taking a simple 
example of the catalytic reaction network with a small number of chemical species. This 
example system is a little specific but can illustrate the changes in the steady distribution 
and the effective network structure consisting only of non-vanishing chemical components. 



(see |2l| for a discreteness- induced switching over states in an autocatalytic network). 

We consider a network of catalytic reactions displayed in Fig.4 (a). Here, we assume that 
ns > always holds. By straightforward calculation, the rate equation for the chemical 
concentrations has a unique fixed-point attractor, which satisfies the relations, 2 < ua >= 
2 < tib >=< nc >=< no >= 2 < tie >= 2 < tir >. 

On the other hand, for small N, some of rii often happens to be 0. In such cases, the 
above relationship on the fixed concentrations no longer holds, where the distribution of the 
chemical concentrations changes temporally among some characteristic distributions with 
different relations among < ua >, < Ub >, < «c >> < n D >, < rig > and < tir >. 

I) State with ua = and n c = 0: When ua happens to be 0, n c also goes to 0, if the 
decrease in nc by the reaction C + E —>■ progresses before the reaction R + C —>■ A + C 
takes place. This state is also reached when n c happens to be and if ha goes to by the 
reaction A + R —> D + R or by the non-catalytic process. Once both ua and rtc vanish, 
this state with extinction of both chemical species is preserved over a long time, because 
neither A nor B is synthesized by catalytic reactions, and only slow non-catalytic changes 
can produce such chemicals. 

At this state, none of the reactions catalyzed by A or C take place. Then, as long as 
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rii > for i ^ A, C, the effective structure of the reaction network is reduced to that 
shown in Fig. 4(b). Only molecules within this sub- network exist. (In Fig. 4(b), the arrows 
indicating the reactions catalyzed by A or C are removed.) In this case, rip increases, 
because the reaction to decrease no catalyzed by A does not take place. Accordingly, the 
reaction S + D — > T + D takes place frequently so that the reaction through the series 
S^T^G^H^ progresses with a high rate. 

II) The state with n G = and n H = 0: When n G happens to be 0, nn also goes to 0, if 
uh decreases by the reaction H + D — >. It also appears when n# happens to be and n G 
decreases to by the reaction G + C — > F + C. Once both n G and n H vanish, this state is 
preserved over a long time. 

Here, when > for % ^ G, H, the effective structure of the reaction network is given in 
Fig. 4(c). Similarly with the case I), chemicals are localized within this subnetwork. In this 
effective network the reaction progresses through the series S-^R^A^, S^R^B^ 
and S — > R — > E — >• with a high rate. 

III) For i y£ A, C, G and H, rii is soon recovered even if it happens to be 0. For example, 
even if rig happens to be 0, as long as no does not reach 0, % is soon recovered by the 
reaction R + D — > B + D. Other chemical species also behave in a similar manner. Then, 
disconnection of the reaction paths does not occur. 

Each of these three states has a long life time when the number of molecules iV is much 
smaller than the number of the chemical species. Stochastic switching over such states is 
commonly observed for a reaction network system with small N, as given in Fig. 2(b). Due 
to these switchings, the temporal fluctuation of each chemical concentration is enhanced, as 
given in Fig. 3(e). Thus, the behavior for small iV is distinct from that obtained in the case 
with iV — > oo. 

Here the relations between <n^>, < n B >, < n G > and < % > differ distinctly from 
that in the continuum limit. For example, < n\ >« tib > and < n G >« no > are 
obtained here. This deviation from the continuum limit is understood easily by considering 
the number distribution at the states I) and II). 

In general, we expect that situations similar to this simple example should appear in some 
part of the reaction networks for a system with a large number of the chemical species, as 
discussed in the last section. This leads to the transition at N < N c , with drastic difference 
in the number distribution from the continuum limit, as given in Fig. 3(a), (b) and (c). 
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5. CRITICAL VALUE OF MOLECULAR NUMBER 



In section 3, we suggested the existence of the critical total number of molecules N c , 
at which the behaviors of random catalytic networks change drastically. For iV > N c , the 
behavior is well represented by continuum description, while for N < N c , deficiency in some 
molecule species suppresses the ongoing reaction. In this section, we study dependence of 
this critical number N c on K and M quantitatively, by computing the reaction rate RR. 

Figure 5 shows examples of dependence of < RR > on N for several values of K, K — 
8, 12, 16, 20 and 24, and M, (a) M = 100 with e = 10~ 4 , (b) M = 1000 with e = 10" 4 . As N 
is decreased, there is a drastic drop in < RR > at some value of N. When N is sufficiently 
larger than this value, < RR > is almost constant, approaching a value at the continuum 
limit ~ 0(< RR >n^oo) where chemical concentrations are almost stationary except for 
small fluctuations. On the other hand, for N smaller than this critical value, the chemical 
concentrations exhibit intermittent switching among several states. 

Now, we consider < RR > for a specific case with N — > oo; for all chemical species, the 
number of the reaction path from other chemical species and those to other chemical species 
are given a unique value K. In such cases, the distribution of the chemical concentrations 
goes to uniform. In each time step, each molecule collides one molecule on average, and 
the probability that the molecule catalyzes the reaction from the collided molecule is K/M. 
Then, < RR > is given by K/M. 

If the fluctuation of the number of the reaction paths for each chemical specie increases, 
the reaction from the chemical species which have small number of reaction paths limit the 
reaction rate. Then, < RR > is smaller than K/M in general. 

Now, to study dependence of < RR > on K, M, and N, we introduce the following 
scaling functions; p =< RR > /(J|) as a normalized reaction rate and NK/M. Note that 
N/M is nothing but the average molecular number of each chemical, which gives the average 
probability that each reaction path is catalyzed. Thus, NK/M gives the average number 
of the effective reaction paths in the network. Based on these considerations, we plot this 
normalized reaction rate as a function of NK/M, to see N, K, M dependence. 

Figure 6 shows p as a function of NK/M for several values of K, K = 8, 12, 16, 20 and 24, 
and M, (a) M = 100, (b) M = 1000 and (c) M = 3000 with e = 10~ 4 . As K is increased, this 
scaled function approaches a form independent of K. Here we note that there is a specific 
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value of NK/M = k c at which all the curves of p for different values of K crosses. Below this 
value k c , the normalized reaction rate decreases with the increase in K, suggesting that the 
deficiency in molecule number per reaction path further suppresses the ongoing reaction. On 
the other hand beyond k c , the normalized reaction rate slightly decreases with the decrease 
in K, due to the increase in the reaction path number fluctuation per path. As the reaction 
is suppressed due to deficiency in molecule number for NK/M < k c , this value of k c gives 
a criterion for discreteness induced transition. 

Dependence of the reaction rate on e is shown in Fig. 7. Here, the value of n c seems to be 
independent of M and e for larger M (M = 1000, M = 3000 or larger M) and smaller e (as is 
also compared with that in Fig.6), where k c ~ 0.8 holds. On the other hand, the value of p at 
NK/M = k c decreases as e — > 0, as shown in Fig. 6(c) and Fig. 7. Then, for NK/M < k c , 
the reaction rate p approaches 0, while p at NK/M > k c is almost unchanged. In other 
words, as e goes to zero, the normalized reaction rate seems to approach a step function with 
for NK/M < k c . This means that the catalytic reactions often freeze. Some reaction paths 
are terminated frequently if NK/M < k c , and only non-catalytic changes give dominant 
contributions to the reaction dynamics. On the other hand, for NK/M > k c , all chemicals 
can react along the connected paths in the catalytic reaction network for most of time, so 
that the behavior at the continuum limit is valid, i.e., < RR >~ 0(< RR >$_ >0O ). To sum 
up the value k c gives a criterion for the discreteness-induced transition, i.e., N c ~ k c x M/K 
for large M. 

6. SUMMARY AND DISCUSSIONS 

In this paper, we have reported a discreteness-induced transition in catalytic reaction 
networks with random connections. When the total number of the molecules is smaller than 
a critical value, transition to a novel dynamical state is observed with a distinct behavior 
from that expected in the continuum limit, i.e., the molecular number — > oo. The behavior 
is characterized by switching over quasi-stationary states where some reaction paths are 
terminated effectively due to deficiency in molecule numbers. Each quasi-stationary state is 
characterized by a smaller set of effective reaction networks. 

The critical molecule number for this transition is shown to be proportional to the number 
of chemical species M divided by the average number of reaction path per species K, with 
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a proportion coefficient, estimated to be k c ~ 0.8 in a limit with large M. Whether the 
number of molecules is large enough to be approximated by the continuous rate equation or 
not is thus determined by N>0.8M/ K. 

So far we have not succeeded in estimating this value of k c analytically. It could possibly 
be related with the percolation threshold of a random network, although the relationship is 
not so straightforward. The percolation transition point k c in general random networks is 
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known to increase logarithmically with the increase in the number of nodes of the network 
while k c decreases to converge to a constant value, with the increase in the number of nodes 
(chemical species). Thus, the relation between k c and k c is still unclear, and analytic estimate 
of k c is still an open question. 

It is noted that the random network does not give a good approximation of the real 
biological network. However, the results obtained in such simplest reaction network should 
give a base to characterize the behaviors in several types of networks. The present study 
gives a starting point for the statistical physics of several catalytic reaction networks. 

On the other hand even by random catalytic reaction networks, some universal features 
such as Zipf 's law in the gene expressions and log normal distributions of the fluctuations of 

n □ 

chemical concentration in cells are reproduced 1 181. [19(. We may expect our results provide 
some insight to biological phenomena. 

In the presented model, below the critical number of molecules (or beyond the critical 
number of molecule species), there appears intermittent transitions over several states, which 
leads to the change in the structure of effective reaction network. Such dynamics may giv e 
a hint to uncover possible mechanisms of switching behavior in the signaling pathways 231 ] . 

Also, dominantly acting paths in the metabolic reaction network of E. coli change in 



accordance with the concentration of the nutrient, as observed by the flux analysis 
In the nutrient poor environment, some chemical species may be deficient in a cell. The 
discreteness induced transition studied here may be relevant to such switch in the effective 
network. 



Recently, 



organisms 
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ar 



ge phenotype fluctuations in isogenetic cells are reported in several 
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261 ] . In such cells, abundances of some molecule species are not so large 



(10 - 100 molecules) while there exist more than 1000 chemical species 271 ] . The intermittent 
switching induced by molecular deficiency may underlie such large fluctuations. 

Last, we briefly discuss dependences of the transition on the network topology or a 
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connecting distribution. We have confirmed that the value of k c as well as the peak (or dip) 
in < rii/N > and ATF are almost unchanged, even if the reaction coefficients are distributed 
uniformly between ands 1, or take two distinct values, say 0.1 and 0.9 randomly. The results 
are not changed either, even if there are some chemical species that catalyze many reaction 
paths and those catalyzing only a few paths. 

However, the value of k c can be shifted if the auto-catalytic productions are dominant. 
If the reaction to a chemical species is catalyzed by itself, its concentration can increase 
rapidly, which may decrease the number of the species at the upstream of this reaction 
towards zero. The probability that some molecule species goes extinct is increased, as has 



been recently reported for a simple reaction network of a few species 

Dependence of the discreteness-induced transition on the network topology is an impor- 
tant future issue both for chemical reaction network dynamics in general and also for the 
understanding of intra-cellular chemical reactions. 
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Awazu-Kaneko Figure 1 



FIG. 1: (a) Illustration of catalytic reaction B + C — > A + C, and (b) an example of catalytic 
reaction networks. 
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Awazu-Kaneko Figure 2 



FIG. 2: Examples of temporal evolutions of for a given reaction network for the case M = 100 
and K = 12 with (a) N = 800 and (b) N = 12. (c) Temporal evolutions of RR for N = 800 
(green) and N = 12 (red). The network is chosen randomly, but the behavior here is typical, and 
observed over most networks generated randomly. 
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FIG. 3: (a)(b) < rtj/iV > of the same system as shown in the previous figure with (a) N = 800 and 
(b) N = 12. (c)(d)(e) N dependency of (c)< rii/N >, (d) the rank of each chemical concentration, 
and (e) < RR > and ATF/(N/M)). 
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FIG. 4: Simple example of catalytic reaction network (a), and effective reaction network consisting 
only of non-vanishing chemicals for the state I) (b) and II) (c) described in the text. Thick and 
dashed arrows indicate the paths with a high reaction rate and those from the extinct chemical 
species at the moment, respectively. 
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FIG. 5: Typical examples of < RR > plotted as a function of N for several K and (a) M = 100 
with e = 10~ 4 and (b) M = 1000 with e = 10~ 4 . 
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FIG. 6: Scaled reaction rate p as a function of NK/M for several values of K and (a) M = 100, 
(b) M = 1000 and (c) M = 3000 with e = 10~ 4 . 
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FIG. 7: Scaled reaction rate p as a, function of NK/M for several values of K with (a) M = 3000 
and e = 10" 5 and (b) M = 3000 and e = 10~ 6 . 
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